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A coarse-grained molecular model, which consists of a spherical particle and an orientation vector, 
is proposed to simulate lipid membrane on a large length scale. The solvent is implicitly represented 
by an effective attractive interaction between particles. A bilayer structure is formed by orientation- 
dependent (tilt and bending) potentials. In this model, the membrane properties (bending rigidity, 
line tension of membrane edge, area compression modulus, lateral diffusion coefficient, and flip-flop 
rate) can be varied over broad ranges. The stability of the bilayer membrane is investigated via 
droplet- vesicle transition. The rupture of the bilayer and worm- like micelle formation can be induced 
by an increase in the spontaneous curvature of the monolayer membrane. 



PACS numbers: 



I. INTRODUCTION 



Amphiphilic molecules, such as lipids and detergents, 
self- assemble into various structures depending on the 
relative size of their hydrophilic parts: spherical or 
worm-like micelles, bilayer membranes, inverted hexago- 
nal structures, and inverted micelles. Among these struc- 
tures, the bilayer membrane of phospholipids has been in- 
tensively investigated, since it is the basic structure of the 
plasma membrane and the intracellular compartments of 
living cells, where the membranes are in a fluid phase 
and lipid molecules can diffuse in quasi-two-dimensional 
space. A vesicle (closed membrane) is considered to be 
a simple model of cells and it has applications in drug- 
delivery systems as a drug carrier. 

Bilayer membranes exhibit many interesting phenom- 
ena such as shape deformation induced by phase sepa- 
ration or chemical reaction, membrane fusion, and mem- 
brane fission. The length scale of these phenomena varies 
from nm to /im, since cells are ~ 10/xm in diameter, 
whereas the thickness of a biomembrane is 5nm. To 
investigate the morphologies of cells and vesicles, the 
molecular structure is assumed to be negligible, and 
the bilayer membrane is described as a smoothly-curved 
mathematical surface^— . The information about the bi- 
layer properties is only reflected in the values of the elas- 
tic parameters. To simulate the membrane with thermal 
fluctuations, a triangulated surface is widely used 4 -! 5 -. An 
alternative model is a meshless membrane^—, where par- 
ticles self-assemble into a membrane by anisotripic poten- 
tial interactions. These models can reproduce /^m-scale 
dynamics of the bilayer membrane well but cannot treat 
a non-bilayer structure such as the stalk structure of a 
membrane fusion intermediat e 14 ' 15 . 

To simulate molecular-scale dynamics and the non- 
bilayer structure, a molecular model is required. Al- 
though computer technology has grown rapidly, the typ- 
ical scale for recent simulations of the all-atom models is 
only 100 ns dynamics of hundreds of lipid molecules. To 
simulate the membranes on longer and larger scales, vari- 
ous coarse-grained molecular models have been proposed 
(see review articles^—) . Recently, the potential param- 



eters in some of the coarse-grained molecular models are 
tuned by atomistic simulations^ 1 - - — . In mapping of inter- 
action parameters, one coarse-grained particle typically 
represents three or four heavy atoms and their accompa- 
nying hydrogen atoms. To further reduce the computa- 
tional costs, larger segments (three or more segment par- 
ticles per amphiphilic molecule) are employed, and the 
solvent is implicitly represented by an effective attrac- 
tive potential between the hydrophobic segment a 19 ' 26 i 2H. 
Model parameters are chosen to generate a bilayer mem- 
brane with reasonably realistic values of elastic proper- 
ties. In this scale, it is difficult to take into account 
chemical details of lipids, such as an unsaturated bond 
in hydrocarbon chains. Instead, this type of models can 
be advantageous to capture the general features in the 
bilayer membrane, since the simplicity of the model can 
allow for wide ranges of variation of the membrane prop- 
erties. 

In this paper, we propose a solvent-free molecular 
model to pursue two purposes: (1) to represent the am- 
phiphilic molecule in a size as small as possible and (2) 
to allow the variation in the membrane properties for 
wide ranges. A molecule consisting of many particles has 
a higher resolution than that with less particles but re- 
quires a smaller length unit and time step for simulations. 
Here, we consider a molecule that consists of a spherical 
particle and an orientation vector. It can reduce compu- 
tational costs to simulate many molecules. In previous 
solvent-free model s 19 ' 25 ~ —, the membrane properties are 
varied only in narrow ranges. On the other hand, in one 
of the meshless membrane models, the bending rigidity 
and the line tension of the membrane edge can be inde- 
pendently varied over wide ranges^. This allows the con- 
ditions of vesicle formation and rupture to be controlled 8 . 
In addition, it is easy to compare the simulation results 
with theoretical predictions. Such tuning capability is 
desired for molecular models. 

In Sec. HH the lipid model and the simulation method 
are described. In Sec. IIII1 the results and discussion 
are provided. The formation of a membrane and its sta- 
bility for droplet- vesicle transitions are described in Sec. 
ITITAl In Sec. HITbI and iHTCl the calculation methods 
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FIG. 1: (Color online) The cutoff function f cu t(rij), the com- 
pact Gaussian weight function w cv (rij), and the repulsive po- 
tential U Tcp (rij). 



and the dependence of the static and dynamic properties 
on model parameters are described, respectively. The 
summary is given in Sec. IIVI 



II. MODEL AND METHOD 

A. Molecular model 

In solvent-free lipid models, an amphiphilic molecule is 
typically represented by three or more particles^.. Here, 
we consider a lipid molecule with minimum (5) degrees 
of freedom for solvent-free molecular simulations. Each 
(i-th) molecule has a spherical particle with an orienta- 
tion vector U; , which represents a direction from the hy- 
drophobic to the hydrophilic part. There are two points 
of interaction in the molecule: the center of a sphere rf 
and a hydrophilic point r° = + u^. The molecules 
interact with each other via the potential, 



U 



^[/ rcp (ry + £ ]Tt/ att (^) (1) 



fetilt 
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where r, 



i"i Ti j — |r^j-|, Ti j — Yi j/vi j, and kj$T 
is the thermal energy. The molecules have an excluded 
volume with a diameter a via the repulsive potential, 



U Isp (r) = exp[-20(r/CT - 1)] 



(2) 



with a cutoff at r — 2. 4a. 

The second term in Eq. (fT|) represents the attractive 
interaction between the molecules. A multibody attrac- 
tive potential U & tt(pi) is employed to mimic the "hy- 
drophobic" interaction. This potential allows the forma- 
tion of the fluid membrane over wide parameter ranges 



and fast lateral diffusion. Similar potentials have been 
applied in the previous membrane models^ ' 26 ' 30 and a 
coarse-grained protein models. The potential U at t{pi) is 
given by 

U att (Pi) = 0-25 ln[l + cxp{-4(p, - p*)}} - C, (3) 

where C = 0.25 ln{l + exp(4p*)} ~ p* is chosen such 
that £7 a tt(0) = 0. The local particle density pi is approx- 
imately the number of particles v\ in the sphere whose 
radius is r att • 



Pi 



/cut(^i j), 



where / cut (r) is a C°° cutoff function, 



fc 



(r) = ( ex P {A(l + (r/rcu ^ ) „_ 1 )} (r < r cut ) 



(r > r cut ) 



(4) 



(5) 



with n = 6, A = ln(2){(r cut /r att ) n - 1}, r att = 1.9a 
(/cutC^att) = 0.5), and the cutoff radius r cut = 2.4ct (see 
Fig. [1]). The potential U at t(pi) acts as a pairwise at- 
tractive potential (U &t t{pi) — P, so that Y,i U &tt{pi) — 
~2 J2i<j /cut(rf j)) for pi < p* - 1 and approaches a con- 
stant value (U a tt(pi) — p*) for pi > p* + 1. It is assumed 
that the hydrophobic parts have no contact with the im- 
plicit solvent (void space) at pi > p* . 



The third and fourth terms in Eq. ([T]) are dis- 
cretized versions of tilt and bending potentials of the tilt 
mode l 32 ' 33 , respectively. A smoothly truncated Gaussian 
function^ is employed as a weight function 



r(r) 



exp( 




(r/Tgaf 

(r/r cc )"- 



j) ( r < r cc) 

(r > r cc ) 



(6) 



with n = 4, r ga = 1.5<7, and r cc = 3a (see Fig. [TJ). All 
orders of derivatives of f C ut(r) and w m \ s (r) are continu- 
ous at the cutoff radii. The weight is a function of rfj 
(not rf ■) to avoid the interaction between the molecules 
in the opposite monolayers of the bilayer. The average 
distance between the neighboring molecules in the same 
monolayer is r n b — 1.05er, and the distance to the neigh- 
boring molecule in the other monolayer is rf^ ~ 3a (see 
Fig. [21a)). Thus, these two potentials act between the 
neighboring molecules in the same monolayer but not 
between the monolayers. The tilt potential has the en- 
ergy minimum in a completely flat membrane with no 
tilt deformation. Similar tilt potentials have been used 
in the meshless membrane models&^r— . The same type 
of bending potential [the fourth term in Eq. |T])] was 
previously used to control the bending rigidity and the 
spontaneous curvature of the monolayer in the molecu- 
lar simulation o i 30 . Positive spontaneous curvature in- 
dicates that the hydrophilic head is larger than the hy- 
drophobic tail of amphiphilic molecules. The bending 
rigidity is numerically calculated and compared with the 
estimation from the continuous description of the mem- 
brane in Sec. IHI CI 
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FIG. 2: (Color online) Probability distribution of (a) the po- 
sitions and (b) the orientation of the molecules in the planar 
membrane at N = 512, p* — 14, e = 2, fcbend — fetilt, and 
Cbd = 0. (a) The z components of rf and rf are shown for 
fetilt = 8. The dashed line represents r| of the molecules at 
u z > 0. (b) The z component u z of the molecular orientation 
is shown for fctilt = 1, 2, 4, and 8. The error bars are displayed 
at several data points. 



B. Brownian dynamics 

We simulated the membrane in the NVT ensemble 
(constant number of molecules N, volume V, and tem- 
perature T) with periodic boundary conditions in a box 
with side length L x , L y , and L z . We employed Brown- 
ian dynamics (molecular dynamics with Langevin ther- 
mostat). The motions of the center of the mass rp = 
(r| + r|)/2 and the orientation Uj are given by under- 
damped Langevin equations: 



drf 
dt 
dv? 



I 



dt 
dijJi 



dt 



~dT = Uh 
-Ccvf +gf(t)+ff, 

-CrU>i + (gj(t) + ff) X + Allj 



(7) 
(8) 
(9) 



where m and / are the mass and the moment of iner- 
tia of the molecule, respectively. The forces are given by 
fp = —dU/drf and f/ = —dU/dui with the perpendic- 
ular component = a — (a • Uj)uj and a Lagrange mul- 
tiplier A to keep uf = 1. According to the fluctuation- 
dissipation theorem, the friction coefficients Cg, Cr an d 
the Gaussian white noises gp(t), (t) obey the follow- 
ing relations: the average (sf^ (£)) = and the variance 



O&Wfl&.W) = 2k B T{ f!l 5 ij 6 aict2 5 01 f )2 6(t-t'), where 
Q ; i,Q;2 € { x iU, z } and Pi, fa € {G,r}. The Langevin 
equations are integrated by the leapfrog algorithm^! with 
Vi,„ = v,(t„) = (vi(t n+1/2 ) +Vi(t n _i/ 2 ))/2. First, the 
velocities are updated by 
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1 i,n+l/2 



a o v ?, n -i/2 + ai(g,> 
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gf (t„)/VAt. 



Then, the positions are updated by 



U i,n+1 
Ui,n+1 



Cr At/27' l + CrAt/2J' 

2^i,n-l/2 ' u i,n 
l+CrAt/2/ ' 



At, 
/2 At, 



= r„. 



i,ra+l/2 



(12) 



J i,n+1 



i.n-\- 



i/K 



n+l I 



We employed to = 1, 7 = 1, Cg = 1, Cr = 1; k B T = 
I, At = 0.005, and the total number of the molecules 
N = 300 to 8192. The results are displayed with the 
length unit a, the energy unit fcsT, and the time unit 
T o — Cg^V^bT 1 - The diffusion coefficient D is normalized 
using the diffusion coefficient Do = a 1 /tq of an isolated 
molecule. The error bars of the data are estimated from 
the standard deviations of three to six independent runs. 



III. RESULTS AND DISCUSSION 
A. Self-assembly and membrane stability 

Molecules self-assemble into spherical droplets at 
^bend = &tilt = 0, i.e., when only the first two terms 
in Eq. ([!} are taken into account. When the third term 
(the tilt potential) is added, the molecules can sponta- 
neously form vesicles. Figure [3] shows the self-assembly 
of molecules from a random gas state. First, small clus- 
ters are formed; these clusters merge into disk-like mi- 
celles. Then, a large disk closes into a vesicle through 
a bowl- like shape (see Fig. EJc)). Similar self-assembly 
processes have been observed in the previous simulations 
of molecular 2 ^ and meshless 8 models. 

In order to clarify the stability of three-dimensional ag- 
gregates and bilayer membranes, the morphologies of the 
aggregates are investigated as fctilt gradually increases or 
decreases. As fctilt increases, a spherical liquid droplet 
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FIG. 3: (Color online) Sequential snapshots of the molecular 
self-assembly at N = 2000, L x = L y = L z = 40a, e = 2, 
p* = 14, fcbcnd = 0, fctnt = 8, and C b d = 0. (a) t/n = 0. (b) 
i/r = 500. (c) t/T = 16650. (d) t/r = 16700. 




FIG. 4: (Color online) Droplet-vesicle transition at N = 500, 
e = 2, p* = 14, fcbcnd = 4, and Cbd = 0. The lower and upper 
lines represent the radius of gyration R g in feyit increasing or 
decreasing, respectively. Sliced snapshots are also shown at 
fctiit = 2.5. 



transforms into a bilayer vesicle. At the transition point 
(fctiit = 4), the radius of gyration R g exhibits an abrupt 
increase as shown in Fig. [4] As fctiit decreases, the 
transition from a vesicle to a droplet occurs, however, 
the transition point (fctiit = 0.5) is much lower. Thus, 
a typical hysteresis for the first-order transition is ob- 
served. The rate of increase or decrease is sufficiently 
low (fctiit = 0.0005t/r ). We checked that the devia- 
tion of the transition points by the annealing rates is 
very small; Afc t iit = 0.1 between knit = 0.000125t/ro 




FIG. 5: (Color online) Formation of vesicles from a droplet 
at N = 4000, e = 2, p* = 14, fc bcn d = 8, and C bd = 0. The 
tilt coefficient ktnt is gradually increased as fc t nt = 0.0005t/ro. 
Sliced snapshots are shown at (a) t/ro = 12500 (fctiit = 6.25), 
(b) t/r = 13600, (c) t/r = 13700, (d) t/r = 14000, and (e) 
t/ro = 14500 (fctiit = 7.25). All molecules are also shown for 
t/r = 13700 in (c'). 



and fc ti i t = 0.001t/ro with N = 500, e = 2, p* = 14, 
fcbcnd = fctiit, and Cbd = 0. The transition points are not 
sensitive to the path of fcbend (fctiit), since the difference of 
the results for fcbcnd = fctiit and constant fcbcnd is smaller 
than their statistical errors. 

For large aggregates with N = 4000, two vesicles or 
a vesicle with disks are formed instead of a single vesi- 
cle. Figure [5] shows an example of the formation of two 
vesicles. A void space is opened in the droplet, and a 
bilayer skirt is formed. Then, it is separated into two 
parts and forms two vesicles. If the separated membrane 
is small, a disk is formed. Since the larger droplets can 
have a clearer molecular layer on the surface, which pre- 
vents the shape change, higher fctiit is needed to trigger 
the shape transition [see Fig. [SJa)]. At N = 300, the co- 
existence region of the droplets and the vesicles is narrow, 
since the number of the molecules is not sufficient to form 
the surface and inside layers. The points of the droplet- 
bilayer transition are also dependent on Cbd, while they 
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FIG. 6: (Color online) Shape transition points of molecu- 
lar aggregates between the droplet and the bilayer membrane 
(vesicles and disks) for various (a) N, (b) Cbd, (c) e, and (d) 
p*. If not specified, N = 500, e — 2, p* = 14, and Cbd = 0. 
The solid and dashed lines represent data with increasing and 
decreasing fctiit, respectively. The open and filled symbols rep- 
resent data with fixed fcbcnd and fcbcnd = fctnt , respectively. 



are almost independent of e and p*. Thus, the bilayer 
stability is determined by the tilt and bending potentials 
but not by the attractive potential. 

Let us discuss the condition required to form a stable 
bilayer. When all molecules have the same orientation 
Uj = Uj, the bending potential energy becomes zero at 
Cbd = 0. Thus, the bending potential with Cbd = 




FIG. 7: (Color online) Surface tension 7 of a flat membrane 
at N = 512, p* = 14, e = 2, fcbcnd = fctnt, and Cbd = 
0. Dependence of 7 on (a) the projected area per molecule 
a xy = 2A xy /Na 2 and (b) the intrinsic area per molecule a = 
2A/Na 2 . The squares, triangles, and circles represent 7 for 
feyit = 2, 4, and 8, respectively. The error bars are smaller 
than the line thickness. 

can have the minimum energy for any structure of the 
aggregate. Therefore, the spherical liquid droplet, which 
has the minimum surface area, would be the equilibrium 
state at fctnt = instead of the bilayer. However, large 
vesicles with N > 1000 and planar membranes can main- 
tain their bilayer structure as a metastable state even at 
fctiit = with finite fcbcnd- Note that the capability to 
keep a pre-formed bilayer membrane does not guarantee 
the self-assembly to the bilayer. In particular, the pe- 
riodic boundary condition is a strong constraint, which 
can keep the bilayer membrane as a thin liquid layer even 
at fctiit = fcbend = 0. In order to obtain the spontaneous 
formation of the bilayer membrane, fctiit > 2 is required. 



B. Calculation of membrane properties 

To investigate the membrane properties, we formed a 
nearly planar membrane without edges or pores. The 
membrane area and the surface tension are varied by in- 
creasing or decreasing the projected area A xy — L x L yi 
where L x = L y . The membrane has a clear bilayer struc- 
ture (see Fig. [5]). The intrinsic area of the tensionless 
membrane per molecule ao — 2Ao/Na 2 , the area com- 
pression modulus Ka, and the half lifetime t$ of the flip- 
flop motion are calculated from the flat membranes with 
TV = 512. The bending rigidity k and the diffusion coef- 
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FIG. 8: (Color online) Spectra of undulation modes {\h(q)\ 2 ) 
of nearly planar, tensionless membranes (7 = 0) at N = 8192, 
p* = 14, e = 2, fcbond = fetilt = 8, and Cbd = 0. Results for 
(\h(q)\ 2 ) calculated from the molecular positions (+) and from 
the averaged positions on a square mesh (x) are shown. The 
inset shows the dependence of 1/ (\h(q)\ 2 ) on q 2 , which is used 
to extract the bending rigidity k. 
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FIG. 9: (Color online) Line tension T of membrane edge at 
e = 2, p" = 14, fcbcnd = 4, fcyit = 4, and C*bd = 0. The circles 
represent Y calculated from a pore on the flat membrane at 
N — 2048. The solid line represents T calculated from the 
striped membrane at N = 512: Fa/k B T = 9.08 ±0.06. 



ficient D are calculated at larger tensionless membranes 
with N = 8192. The line tension V of membrane edge 
is calculated from the strip of the flat membrane with 
N = 512. 

The surface tension 7 is given h y 34 i 35 

7 - (Pzz - (P XX + Pyy)/2)L Z , (13) 

with the diagonal components of the pressure tensor 

, dU . 



P aa = {Nk B T 



(14) 



where a 6 {x, y, z}. When the potential interaction 
crosses the periodic boundary, the periodic image ai + 
nh a nearest to the other interacting molecules is em- 
ployed. The intrinsic area A of the membrane is larger 
than the projected area A xy in the xy-plane due to the 
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FIG. 10: (Color online) Time development of the probability 
difference P up — -Pdown of the molecules in the upper and lower 
monolayers at TV = 512, p* — 14, e — 2, febend — 0, and 
Cbd = 0. At the initial states (t = 0), P up = 1 and Pdown = 0. 
The error bars are displayed at several data points. 



membrane undulations. We calculate A from a y/N/2 x 

\/N/2 square mesh with (x mb ,y mh ) = (d m h«, dmhj)- The 
height z m h of a mesh point is obtained from the weighted 
average of molecular position in the four neighbor 
cells, with z mh = £\ z i w nlh {x l , y t )/(J2i w^h{x l7 y t )) and 

W m b(Xi,yi) = (1 - \Xi - X mh |/rf mh )(l - \Vi - 2/mh|/rfmh)- 

Figure [7] shows the area dependence on the surface ten- 
sion 7. The tension 7 exhibits a roughly linear increase 
with the molecular area at 7 > 0. The compressed mem- 
brane with 7 < buckles out of plane and has the larger 
intrinsic area A than the projected area A xy . Similar 7 
dependence and buckling are obtained in the simulations 
of other molecular model s 29 ' 36 and meshless models^^. 

The area Ao of the tensionless membrane (7 = 0) is 
obtained by the minimization of 7, where the projected 
area is updated as A xy w = A xy — i"yAi 7 every At 7 in- 
terval, where 7 is the time average for Af 7 . We use 
6t = 0.00025 to 0.005 and Ai 7 /V = 50 or 100. The 
area compression modulus is defined as 



K A = A dj/dA\ A=Ao . 



(15) 



We calculate Ka from the slope of a-7 lines shown in Fig. 
db). 

The bending rigidity k is calculated from the spectra 
of undulation modes (|ft-((?)| 2 ) of the planar membranes 
in Fourier space d 37 ' 38 , 



(\Hq)\ 2 ) = , 
79 



(16) 



Figure [8] clearly shows the q~ A dependence of the ten- 
sionless membrane. We calculate |/i(<?)| 2 from the raw 
data (the particle position r|), as well as from the square 
mesh with the same mesh-points which were used for the 
estimation of the intrinsic area A. Averaging over the 
mesh removes most of the effects of the molecular pro- 
trusions. The bending rigidity k is estimated from a fit 
of l/(\h(q)\ 2 ) = {K/k B T)(q 2 ) 2 for (q/n) 2 < 0.015, where 
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the difference of two spectra is very small (see the inset 
of Fig.©. 

The line tension T of the membrane edge is calculated 
from the strip of the flat membrane a a 39 ' 40 

T = (P xx - (Pyy + P zz )/2)L y L z /2, (17) 

since T is the energy per unit length of the membrane 
edge, and the length of the membrane edge is 2L X . Since 
the striped membrane is tensionless, (P yy ) — (Pzz) — 
for solvent-free simulations. The tension T and its error 
bar are estimated from the average and standard devi- 
ations for L x /a = 14, 15, 16, and 18. Alternatively, T 
can be also calculated from a circular pore on the flat 
membraneii^ 9 .. In this case, T is balanced with the sur- 
face tension 7 as T — 7r p o. Since one has to estimate the 
pore radius r p0 , this method gives larger statistical errors 
as shown in Fig. |H1 Therefore, we used the membrane 
strip for the calculation of T. 

In bilayer membranes, molecules can move laterally on 
a monolayer and transversely between upper and lower 
monolayers. The lateral diffusion coefficient D of the 
molecules is calculated from the diffusion of the molecu- 
lar projections in the xy plane; D — ((xi(t) — Xi(0)) 2 + 
(Vi(t) ~ yi(0)) 2 )/4t. In the fluid phase, the molecules 
exhibit a fast diffusion rate D/Do = 0.05 to 0.1. 

The relaxation time of the transverse motion (flip-flop) 
between the upper and lower monolayers is measured 
from the relaxation of the labeled molecule s 26 ' 41 . The dif- 
ferential equation of the probability P np (t) (Pdown (i)) of 
the molecules, which belong to the upper (lower) mono- 
layer, is given by dP up /dt = ~k u P up + k d P down , where 
■Pup + Pdown = 1- For planar membranes, fc u = fed, and 
P up = Pdown = 1/2 at i — >oo. When the molecules in 
the upper monolayer are initially labeled (P up (0) = 1), 
the probability decays as 

PupW - Pdow„(t) = exp[-(fc u + k d )t] (18) 

with the half lifetime = ln(2)/(fc u + fed). Figure [TU] 
shows that P up (t) — Pdown (t) indeed follows the exponen- 
tial decay in our simulations. Either the z component of 
the position rf or orientation can be used to detect the 
nomolayer, to which a molecule belongs. Both r° and 
give the same probability distribution (entirely same for 
most of the parameters) , since both the quantities have 
a clear minimum between two peaks (see Fig. [2]). 

C. Parameter dependence of membrane properties 

Figures ITTHT51 show the parameter dependence of the 
properties of the tensionless membrane. The bilayer 
membrane is formed in the fluid phase over broad ranges 
of the parameters due to the multibody attractive poten- 
tial. A gel phase is obtained only at a large value of p* in 
Eq. ([3]). The fluid-gel transition occurs at p* = 16 and 
p* = 15 for e = 2 and 8, respectively [see jumps of D in 
Fig. UTT e)]. As p* decreases, the lateral diffusion and flip- 
flop motion become faster, and the membrane elasticities 




12 14 „, 16 



FIG. 11: (Color online) Parameter p* dependence of (a) the 
intrinsic area 00 = 2Ao/Na 2 per molecule, (b) area compres- 
sion modulus Ka, (c) bending rigidity k, (d) line tension F, 
and (e) diffusion coefficient D for the tensionless membrane 
at fctiit = 4, fcbend = 4, and Cbd = 0. The circles and squares 
represent data for e = 2 and 8, respectively. 



(Ka, k, T) decrease [see Figs. HTIandfTSFd)]. The intrin- 
sic area an. = 2Ao/Na 2 per molecule (N/2 molecules in 
one of monolayers) decreases with increasing p* and ap- 
proaches the closest-packing area v / 3cr 2 /2 ~ 0.87er 2 in 
the gel phase. In this paper, we focus on the fluid mem- 
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FIG. 12: (Color online) Parameter e dependence of (a) 
a = 2A /Na 2 , (b) K A , (c) K, (d) Y, and (e) D for p* = 14, 
fcbond = fctut, and Cbd = 0. The triangles, circles, and squares 
represent data for fctut = 2, 4, and 8, respectively. 



brane and set p* — 14, hereafter. 

The dependence on the strength of attraction e in Eq. 
([1]) is shown in Fig. [121 It has a tendency similar to p* 
dependence. The line tension Y can be varied by e. At 
e = 1, r is close to k^T/a and the bilayer membrane 
is accompanied by free molecules (gas) with the average 
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FIG. 13: (Color online) Bending rigidity k dependence on (a) 
fctut, (b) fcbend, and (c) fctut + fcbend for p* = 14, e = 2, and 
Cbd = 0. The solid lines with squares, circles, and triangles 
represent data for fc bcnd = 0, fcbend = fctut, and fc bend = 8, 
respectively. The dashed lines with crosses and diamonds 
represent data for fctut = 4 and 8, respectively. 



density of the gas ~ 0.001/ct 3 . The molecules depart 
from the bilayer membrane and return. At e = 0.75 with 
fcbend = fctnt = 4 and Cbd = 0, the bilayer membrane 
breaks and small micelles are formed (~ 60 molecules 
per micelle for N/V = 0.08/er 3 ). On the other hand, no 
free molecules are seen at e > 2. The critical micelle con- 
centration (CMC) of lipids is very low, and their chemical 
potential difference in solution and in membrane is typ- 
ically more than lOfce? 1 per lipid*^. Thus, the number 
of lipid molecules on vesicles is conserved in typical ex- 
periments. In order to keep the number of molecules on 
membrane constant during simulations, a sufficiently low 
CMC is required. We mainly use e = 2 and 8, where the 
fluid membranes without free molecules are obtained. 

The bending rigidity k of the bilayer membrane can 
be controlled by fctnt or fcbend- in Eq. ([T]). Figures 
I13f a) and (b) show the linear dependence of k on fc t nt 
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FIG. 14: Parameter fctat dependence of (a) ao = 2Ao/Na 2 , 
(b) K A , (c) r, and (d) D for p* = 14, s = 2, and C bd = 0. The 
squares, circles, and triangles represent data for febend = 0, 
febend = fetut, and febend = 8, respectively. 
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FIG. 15: (Color online) Parameter febend dependence of (a) 
a = 2A /Na 2 , (b) K A , (c) T, and (d) D for p* = 14, e = 2, 
and Cbd = 0. The squares and circles represent data for 
fetut = 4 and 8, respectively. 



and febend, respectively. When they are plotted together 
for fctiit + fcbcnd, ah lines are roughly overlapped on a 
line K/k B T = 2(fctat + febend) + 6 [see Fig. Eftc)]. A 
large deviation from the line is seen only at one data 
point at fetut = and febend = 8 (the leftmost triangle), 
where the bilayer structure is metastable. The bend- 
ing rigidity k is also weakly dependent on e. The k-e 
curve maintains its shape and shifts upward with in- 
creasing fctiit as shown in Fig. I12f c). Thus, it is ex- 
pressed by n/k B T = 2(fc t iit + &bcnd) + b s (e), where b £ (e) 
is an increasing function as b £ (2) — 6, 6 e (4) = 14, and 
6 e (8) = 18. The area compression modulus Ka increases 
with increasing fetut j while Ka shows only slight depen- 
dence on febend for large fetut (see Figs. ITU and [To|) . The 
other membrane properties a, T, and D show weak de- 



pendence on fetut and febend- Thus, k can be varied by 
fetut or febend without a large variation in T. The modulus 
Ka can be varied by fetut- 

The bending elastisity generated by the orientation- 
dependent potentials can be derived from the Helfrich 
theory for monolayer membranes. When the orientation 
vectors are equal to the normal vectors of the monolay- 
ers without tilt deformation, the bending and tilt energies 
are written by 

Ucv = JdA ^[(<7i - Q 2 + (C a - C ) 2 } 

+%(C?+Cf)] (19) 
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FIG. 16: (Color online) Line tension F dependence on Cbd 
for (a) fcbend = fctut and (b) fcbcnd = 8 at p* — 14 and e — 2. 
(a) The diamonds, squares, triangles, and circles represent 
data for fc t nt = 2, 4, 6, and 8, respectively, (b) The squares, 
triangles, and circles represent data for fctut = 1,2, and 8, 
respectively. 

= J dA "bend + <m [{Ci +C2 _ Co) 2 

-« C n d + <nt)CiC 2 + U (20) 

in the continuum limit, where C\ and Ci are two prin- 
cipal curvatures of the monolayer. The first and second 
terms in Eq. (|T5|) are the contributions of the bending 
and tilt potentials, respectively. The spontaneous curva- 
ture of the bending potential is given by C — Cbd/^nb^ 
The nearest-neighbor distance f n b — 1.05u is obtained 
from the radial distribution function. By assuming the 
hexagonal packing of the molecules in the monolayers, 
the monolayer bending rigidities generated by the bend- 
ing and tilt potentials are estimated as K^^/k^T = 
bcnd^cv (r nb ) and n' tnt /k B T = V3k t g t w cv (r nh )/2, re- 
spectively. The bending rigidity of the monolayer is given 
by the sum of these K mono = K bond + n' tilt . For the 
monolayer membrane, Eq. (|2TJ)) gives the saddle-splay 
modulus Smono = -Kmono and the spontaneous curva- 
ture C = {^bcnd/Kcnd + <ilt)} C 'bd/r„b with U = 

Kend + «4it)(V 2 + <iit/ K bend) c o- Thus , the bending 
rigidity k of the bilayer is estimated as k — 2/c mono ~ 
(2.1fc b0 „d + l.lfetiit)fcBT from w cv (1.05er) = 0.61. This 
estimation of K supports the linear dependence of the 
obtained simulation results on fcbcnd and fctut- The pref- 
actor (~ 2) of fcbend gives the quantitative agreement, 
whereas the prefactor of fc t u t is half of the numerical es- 
timation. In the simulation, the thermal fluctuations in- 




FIG. 17: (Color online) Parameter Cbd dependence of (a) 
do = 2A /Na 2 , (b) K A , (c) k, and (d) D for p* = 14, e = 2, 
and fcbond = fctiit- The squares and circles represent data for 
fctiit = 4 and 8, respectively. 



duce molecular protrusion and tilt. These tilt fluctua- 
tions likely change the prefactor of fctut to twice its value 
(Kj ilt /K bend = fctut /fcbend)- The attractive potential also 
adds a small bending resistance, b e {e). 

The flip-flop time Tg shows exponential dependence on 
the parameters of the tilt and bending potentials, while it 
has weak dependence on e (see Fig. IT51) . The free-energy 
barrier between two monolayers would be the main fac- 
tor to determine Tg. It can be roughly estimated from 
the probability distribution of the molecular orientation 
as Fg = k B T(\n[P(u z = 1)] - \n[P(u z = 0)]) [see Fig. 
[2jb)]. The dashed lines in Fig. [LSta) represent Fg cal- 
culated by this method. It gives very good agreements 
with In(rgt/ro), and thus, the flip-flop time is written by 
T ff/ T o = bg exp(Fg / k B T) with bg ~ 1. The barrier Fg is 



11 




1 1 1 1 


1 | 1 1 1 | 


1 1 
-® 4 





jgf^ k tilt = 4 
>^ k bend = 4 





r i . . 


i 1 i i i 1 


(b) - 


-0.4 


o r 0.4 






10" - 



10° — 



1(T — 



-® # 

k tilt = 4 



k tilt - 2 



i ■ ■ ■ r 
g — ®- 



(C) E 



e B B 10 




FIG. 18: (Color online) Half lifetime tb of flip-fiop motion. 

(a) Dependence on fc t ilt at p* = 14, e — 2, and Cbd = 0. 
The dashed lines represent the free-energy barrier Fg/ksT 
estimated by the orientation distribution shown in Fig. [2jb). 

(b) Dependence on Cbd at p* = 14 and e = 2. (c) Dependence 
on e at p* — 14, fcbcnd = fctat> an d Cbd = 0. (d) Dependence 
on p* at fcbcnd = fctiit = 4 and C b d = 0. 




FIG. 19: (Color online) Sequential snapshots of vesicle rup- 
ture at AT = 2000, e = 2, p* = 14, fcbend = 8, fctiit = 8, and 
Cbd = 0.85. (a) t/ro = 500. (b) t/r = 640. (c) t/r = 3000. 



typical simulation time scales. This result agrees with 
the experimental observations. On the other hand, at 
small fctiit and fcbcnd, the molecule shows very fast flip- 
flop, which is advantageous to equilibrate the membrane 
system quickly during simulations. In the present model, 
one can choose slow or fast flip-flop condition to match 
one's simulation purpose. 

The spontaneous curvature of the monolayer is var- 
ied by the parameter Cbd- According to the above con- 
tinuous theory, C = {k hcnd /a(k hcnd + fctiit)}C b d- At 
high spontaneous curvature Cq ~ 1/er, a worm-like mi- 
celle is stabilized. As Cbd increases, the line tension T 
decreases, and the bilayer structure becomes unstable 
at T ~ k^T/a. The Cbd~r curves in Fig. [TBI have a 
parabolic shape with a maximum at Cbd ~ —0.1. The 
bending rigidity k and the flip-flop time increase with 
Cbd, while ao, Ka, and D do not depend significantly 
on Cbd [see Figs. H71 and H5T b)]. At negative Co, the 
molecules easily stay at the middle of the bilayer struc- 
ture (low Fff), and the less clear bilayer would generate 
lower k. 

Membrane rupture is observed when Cbd is increased. 
Figure [TH] shows the rupture process of a vesicle at Cbd = 
0.85. The initial vesicle is metastable and maintains its 
shape for 500to. Then, a pore opens and grows to cracks 
[see Fig. HW b)]. Finally, a branched worm-like micelle is 
formed [see Fig. [TWc)]. Thus, the line tension and the 
stable structure can be varied by Cbd- 



linear with the orientation-dependent potentials, while it 
is almost independent of the attractive potentials. The 
flip-flop rate can be tuned by fctiit or fcbcnd- In typical 
experimental conditions, the flip-flop motion of phospho- 
lipids is very slow, and rg is several hours or days 41 . 
At a sufficiently high k t nt or fcbcnd, no flip-flop occurs in 



IV. SUMMARY 

Wc have proposed a simple molecular model of bilayer 
membranes. The molecule has five degrees of freedom, 
three translational degrees and two orientational degrees. 
Since this molecular model consists of one spherical par- 
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ticle for the excluded volume, it is smaller than previous 
molecular models. Thus, this model is more efficient for 
large-scale simulations. Despite the extreme simplifica- 
tion, this model reproduces many aspects of lipid bilayer 
membranes. 

In the present model, the properties of the fluid 
membranes can be controlled in broad ranges including 
metastable bilayer membranes. The bending rigidity k is 
linearly dependent on k t nt and fcb cn d- The line tension V 
of the membrane edge can be varied by e and Cbd ■ The 
membrane has a wide range of fluid phase, and the fluid- 
gel transition point can be controlled by p*. The area 
compression modulus if a can be varied by fctnt- The 
flip-flop time Tg can be varied by fctilt and fcbend- 

The corresponding time and length scales can be 
mapped by the membrane thickness 5nm and the lateral 
diffusion coefficient ~ 10~ 8 cm 2 /s for phospholipids 43 . 
Thus, the unit length and times are estimated as a = 2 
nm and To ~ 0.1/^s, respectively. 

Our model is suitable to study the details of 
topological-change processes of the membrane, such as 
molecular self-assembly, pore formation, membrane fu- 
sion, and membrane fission. Although meshless mem- 
brane models can also be used, they cannot treat detailed 
structures such as the fusion intermediates 1 - 4 * 1 ^. Re- 
cently, fusion dynamics have been investigated by molec- 
ular simulations 4 ^ - — . However, the condition to deter- 
mine the fusion pathways has not be clarified so far. 



The ability to vary the membrane properties in wide 
ranges would be an advantage of this model for quan- 
titative investigation of the membrane fusion pathways. 
On the other hand, this model is not suitable for ap- 
plication to phenomena, in which the stretching of hy- 
drophobic chains or atomistic details play an important 
role, e.g., in the interactions of membrane proteins via 
hydrophibic mismatch, where molecular stretching is not 
negligible^ - — . 

In this paper, we used Brownian dynamics but one 
can also use the Monte Carlo method and molecular dy- 
namics with another thermostat. When the solvent- free 
molecular model is combined with a particle-based hy- 
drodynamics method, multi-particle collision dynamics 
(MFC) 56 ' 5 , the hydrodynamic interaction can be taken 
into account as demonstrated for a meshless membrane 
models. Thus, the present coarse-grained molecular 
model is efficient for large-scale dynamics of a biomem- 
brane with or without hydrodynamic interactions and is 
applicable to many kinds of phenomena. 
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